{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to atomman: ElasticConstants class\n",
    "\n",
    "__Lucas M. Hale__, [lucas.hale@nist.gov](mailto:lucas.hale@nist.gov?Subject=ipr-demo), _Materials Science and Engineering Division, NIST_.\n",
    "    \n",
    "[Disclaimers](http://www.nist.gov/public_affairs/disclaimer.cfm) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction<a id='section1'></a>\n",
    "\n",
    "The ElasticConstants class represents the elastic constants of a crystal. The class methods focus on:\n",
    "\n",
    "- Allowing values to be set/retrieved in a number of different formats.\n",
    "- Correctly handling transformations to other Cartesian orientations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Library Imports**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "atomman version = 1.4.0\n",
      "Notebook executed on 2021-08-05\n"
     ]
    }
   ],
   "source": [
    "# Standard Python libraries\n",
    "import datetime\n",
    "\n",
    "# http://www.numpy.org/\n",
    "import numpy as np\n",
    "np.set_printoptions(precision=4, suppress=True)\n",
    "\n",
    "# https://github.com/usnistgov/atomman\n",
    "import atomman as am\n",
    "import atomman.unitconvert as uc\n",
    "\n",
    "# Show atomman version\n",
    "print('atomman version =', am.__version__)\n",
    "\n",
    "# Show date of Notebook execution\n",
    "print('Notebook executed on', datetime.date.today())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Elastic constants representations<a id='section2'></a>\n",
    "\n",
    "The ElasticConstants object allows for various representations of the elastic constants to be retrieved:\n",
    "\n",
    "- __Cij__ (6, 6) array of Voigt representation of elastic stiffness.\n",
    "- __Sij__ (6, 6) array of Voigt representation of elastic compliance.\n",
    "- __Cij9__ (9, 9) array representation of elastic stiffness.\n",
    "- __Cijkl__ (3, 3, 3, 3) array representation of elastic stiffness.\n",
    "- __Sijkl__ (3, 3, 3, 3) array representation of elastic compliance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1. Build demonstration (see [Section #3](#section3) for more setting options)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 0 K values for Al\n",
    "C11 = uc.set_in_units(1.143e12, 'dyn/cm^2')\n",
    "C12 = uc.set_in_units(0.619e12, 'dyn/cm^2')\n",
    "C44 = uc.set_in_units(0.316e12, 'dyn/cm^2')\n",
    "\n",
    "C = am.ElasticConstants(C11=C11, C12=C12, C44=C44)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2. Show different representations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cij** is the 6x6 Voigt representation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cij (GPa) ->\n",
      "[[114.3  61.9  61.9   0.    0.    0. ]\n",
      " [ 61.9 114.3  61.9   0.    0.    0. ]\n",
      " [ 61.9  61.9 114.3   0.    0.    0. ]\n",
      " [  0.    0.    0.   31.6   0.    0. ]\n",
      " [  0.    0.    0.    0.   31.6   0. ]\n",
      " [  0.    0.    0.    0.    0.   31.6]]\n"
     ]
    }
   ],
   "source": [
    "print('Cij (GPa) ->')\n",
    "print(uc.get_in_units(C.Cij, 'GPa'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cij9** is the full 9x9 representation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cij9 (GPa) ->\n",
      "[[114.3  61.9  61.9   0.    0.    0.    0.    0.    0. ]\n",
      " [ 61.9 114.3  61.9   0.    0.    0.    0.    0.    0. ]\n",
      " [ 61.9  61.9 114.3   0.    0.    0.    0.    0.    0. ]\n",
      " [  0.    0.    0.   31.6   0.    0.   31.6   0.    0. ]\n",
      " [  0.    0.    0.    0.   31.6   0.    0.   31.6   0. ]\n",
      " [  0.    0.    0.    0.    0.   31.6   0.    0.   31.6]\n",
      " [  0.    0.    0.   31.6   0.    0.   31.6   0.    0. ]\n",
      " [  0.    0.    0.    0.   31.6   0.    0.   31.6   0. ]\n",
      " [  0.    0.    0.    0.    0.   31.6   0.    0.   31.6]]\n"
     ]
    }
   ],
   "source": [
    "print('Cij9 (GPa) ->')\n",
    "print(uc.get_in_units(C.Cij9, 'GPa'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cijkl** is the full 3x3x3x3 representation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cijkl (GPa) ->\n",
      "[[[[114.3   0.    0. ]\n",
      "   [  0.   61.9   0. ]\n",
      "   [  0.    0.   61.9]]\n",
      "\n",
      "  [[  0.   31.6   0. ]\n",
      "   [ 31.6   0.    0. ]\n",
      "   [  0.    0.    0. ]]\n",
      "\n",
      "  [[  0.    0.   31.6]\n",
      "   [  0.    0.    0. ]\n",
      "   [ 31.6   0.    0. ]]]\n",
      "\n",
      "\n",
      " [[[  0.   31.6   0. ]\n",
      "   [ 31.6   0.    0. ]\n",
      "   [  0.    0.    0. ]]\n",
      "\n",
      "  [[ 61.9   0.    0. ]\n",
      "   [  0.  114.3   0. ]\n",
      "   [  0.    0.   61.9]]\n",
      "\n",
      "  [[  0.    0.    0. ]\n",
      "   [  0.    0.   31.6]\n",
      "   [  0.   31.6   0. ]]]\n",
      "\n",
      "\n",
      " [[[  0.    0.   31.6]\n",
      "   [  0.    0.    0. ]\n",
      "   [ 31.6   0.    0. ]]\n",
      "\n",
      "  [[  0.    0.    0. ]\n",
      "   [  0.    0.   31.6]\n",
      "   [  0.   31.6   0. ]]\n",
      "\n",
      "  [[ 61.9   0.    0. ]\n",
      "   [  0.   61.9   0. ]\n",
      "   [  0.    0.  114.3]]]]\n"
     ]
    }
   ],
   "source": [
    "print('Cijkl (GPa) ->')\n",
    "print(uc.get_in_units(C.Cijkl, 'GPa'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Sij** is the Voigt 6x6 elastic compliances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sij (1/GPa) ->\n",
      "[[ 0.0141 -0.005  -0.005   0.      0.      0.    ]\n",
      " [-0.005   0.0141 -0.005   0.      0.      0.    ]\n",
      " [-0.005  -0.005   0.0141  0.      0.      0.    ]\n",
      " [ 0.     -0.     -0.      0.0316 -0.     -0.    ]\n",
      " [ 0.      0.      0.      0.      0.0316  0.    ]\n",
      " [ 0.      0.      0.      0.      0.      0.0316]]\n"
     ]
    }
   ],
   "source": [
    "print('Sij (1/GPa) ->')\n",
    "print(uc.get_in_units(C.Sij, '1/GPa'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Sijkl** is the full 3x3x3x3 elastic compliances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sijkl (1/GPa) ->\n",
      "[[[[ 0.0141  0.      0.    ]\n",
      "   [ 0.     -0.005   0.    ]\n",
      "   [ 0.      0.     -0.005 ]]\n",
      "\n",
      "  [[ 0.      0.0079  0.    ]\n",
      "   [ 0.0079  0.      0.    ]\n",
      "   [ 0.      0.      0.    ]]\n",
      "\n",
      "  [[ 0.      0.      0.0079]\n",
      "   [ 0.      0.      0.    ]\n",
      "   [ 0.0079  0.      0.    ]]]\n",
      "\n",
      "\n",
      " [[[ 0.      0.0079  0.    ]\n",
      "   [ 0.0079  0.      0.    ]\n",
      "   [ 0.      0.      0.    ]]\n",
      "\n",
      "  [[-0.005   0.      0.    ]\n",
      "   [ 0.      0.0141  0.    ]\n",
      "   [ 0.      0.     -0.005 ]]\n",
      "\n",
      "  [[ 0.     -0.     -0.    ]\n",
      "   [-0.     -0.      0.0079]\n",
      "   [-0.      0.0079 -0.    ]]]\n",
      "\n",
      "\n",
      " [[[ 0.      0.      0.0079]\n",
      "   [ 0.      0.      0.    ]\n",
      "   [ 0.0079  0.      0.    ]]\n",
      "\n",
      "  [[ 0.     -0.     -0.    ]\n",
      "   [-0.     -0.      0.0079]\n",
      "   [-0.      0.0079 -0.    ]]\n",
      "\n",
      "  [[-0.005   0.      0.    ]\n",
      "   [ 0.     -0.005   0.    ]\n",
      "   [ 0.      0.      0.0141]]]]\n"
     ]
    }
   ],
   "source": [
    "print('Sijkl (1/GPa) ->')\n",
    "print(uc.get_in_units(C.Sijkl, '1/GPa'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Setting values<a id='section3'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1. Setting using attributes\n",
    "\n",
    "The values of the elastic constants can be initialized or set using any of the above representations.\n",
    "\n",
    "- During initialization, pass one of the following as the only parameter: Cij, Sij, Cij9, Cijkl, or Sijkl.\n",
    "- For an already initialized object, set any of the representations directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "Sijkl = C.Sijkl\n",
    "\n",
    "# Initialize a new C with Sijkl values\n",
    "newC = am.ElasticConstants(Sijkl = Sijkl)\n",
    "\n",
    "# Show values to be the same\n",
    "print(np.allclose(newC.Cij, C.Cij))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0.]]\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "# Initialize all zeros elastic constants array\n",
    "newC = am.ElasticConstants()\n",
    "print(newC)\n",
    "\n",
    "# Set values using Cij9\n",
    "newC.Cij9 = C.Cij9\n",
    "\n",
    "# Show values to be the same\n",
    "print(np.allclose(newC.Cij, C.Cij))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2. Setting using crystal-specific constants\n",
    "\n",
    "Alternatively, the elastic constants can be defined by supplying a full set of unique constants for a given crystal structure in the standard reference frame. During initialization, the specific form is inferred from the given parameters. After initialization, the functions corresponding to a specific crystal family can be called. \n",
    "\n",
    "- **isotropic** (two unique values required)\n",
    "    - **C11**\n",
    "    - **C12** \n",
    "    - **C44** (2\\*C44 = C11-C12)\n",
    "    - **M** P-wave modulus (M = C11)\n",
    "    - **lambda** Lame's first parameter (lambda = C12)\n",
    "    - **mu** Shear modulus (mu = C44)\n",
    "    - **E** Young's modulus\n",
    "    - **nu** Poisson's ratio\n",
    "    - **K** Bulk modulus\n",
    "- **cubic** (all three values required)\n",
    "    - **C11**\n",
    "    - **C12**\n",
    "    - **C44**\n",
    "- **hexagonal** (five unique values required) \n",
    "    - **C11**\n",
    "    - **C12**\n",
    "    - **C13**\n",
    "    - **C33**\n",
    "    - **C44**\n",
    "    - **C66** (2\\*C66 = C11-C12)\n",
    "- **tetragonal** (six values required, one optional)\n",
    "    - **C11**\n",
    "    - **C12**\n",
    "    - **C13**\n",
    "    - **C16** optional (C16=0 for some space groups)\n",
    "    - **C33**\n",
    "    - **C44**\n",
    "    - **C66**\n",
    "- **rhombohedral** (six unique values required, one optional)\n",
    "    - **C11**\n",
    "    - **C12**\n",
    "    - **C13**\n",
    "    - **C14**\n",
    "    - **C15** optional (C15=0 for some space groups)\n",
    "    - **C33**\n",
    "    - **C44**\n",
    "    - **C66** (2\\*C66 = C11-C12)\n",
    "- **orthorhombic** (all nine values required)\n",
    "    - **C11**\n",
    "    - **C12**\n",
    "    - **C13**\n",
    "    - **C22**\n",
    "    - **C23**\n",
    "    - **C33**\n",
    "    - **C44**\n",
    "    - **C55**\n",
    "    - **C66**\n",
    "- **monoclinic** (all thirteen values required)\n",
    "    - **C11**\n",
    "    - **C12**\n",
    "    - **C13**\n",
    "    - **C15**\n",
    "    - **C22**\n",
    "    - **C23**\n",
    "    - **C25**\n",
    "    - **C33**\n",
    "    - **C35**\n",
    "    - **C44**\n",
    "    - **C46**\n",
    "    - **C55**\n",
    "    - **C66**\n",
    "- **triclinic** (all twenty-one values required)\n",
    "    - **all Cij where i <= j**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cij (GPa) ->\n",
      "[[144.  74.  67.   0.   0.   0.]\n",
      " [ 74. 144.  67.   0.   0.   0.]\n",
      " [ 67.  67. 166.   0.   0.   0.]\n",
      " [  0.   0.   0.  33.   0.   0.]\n",
      " [  0.   0.   0.   0.  33.   0.]\n",
      " [  0.   0.   0.   0.   0.  35.]]\n",
      "\n",
      "Does 2*C66 = C11-C12?\n",
      "C11 - C12 (GPa) = 70.0\n",
      "2 * C66 (GPa) =   70.0\n"
     ]
    }
   ],
   "source": [
    "# Define dict with Zr HCP constants\n",
    "Cdict = {}\n",
    "Cdict['C11'] = uc.set_in_units(144.0, 'GPa')\n",
    "Cdict['C12'] = uc.set_in_units(74.0, 'GPa')\n",
    "Cdict['C13'] = uc.set_in_units(67.0, 'GPa')\n",
    "Cdict['C33'] = uc.set_in_units(166.0, 'GPa')\n",
    "Cdict['C44'] = uc.set_in_units(33.0, 'GPa')\n",
    "\n",
    "# Initialize by passing dict key-values as parameters\n",
    "C = am.ElasticConstants(**Cdict)\n",
    "\n",
    "# Show that Cij array is properly constructed\n",
    "print('Cij (GPa) ->')\n",
    "print(uc.get_in_units(C.Cij, 'GPa'))\n",
    "print()\n",
    "\n",
    "# Show that 2 * C66 = C11 - C12\n",
    "print('Does 2*C66 = C11-C12?')\n",
    "print('C11 - C12 (GPa) =', uc.get_in_units(C.Cij[0,0] - C.Cij[0,1], 'GPa'))\n",
    "print('2 * C66 (GPa) =  ', uc.get_in_units(2 * C.Cij[5,5], 'GPa'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cij (GPa) ->\n",
      "[[ 74.36   26.19   29.17   15.41    0.      0.   ]\n",
      " [ 26.19   74.36   29.17  -15.41    0.      0.   ]\n",
      " [ 29.17   29.17   51.6     0.      0.      0.   ]\n",
      " [ 15.41  -15.41    0.     31.35    0.      0.   ]\n",
      " [  0.      0.      0.      0.     31.35   15.41 ]\n",
      " [  0.      0.      0.      0.     15.41   24.085]]\n",
      "\n",
      "Does 2*C66 = C11-C12?\n",
      "C11 - C12 (GPa) = 48.170000000000016\n",
      "2 * C66 (GPa) =   48.170000000000016\n"
     ]
    }
   ],
   "source": [
    "# Define dict with Bi2Te3 rhombohedral constants\n",
    "Cdict = {}\n",
    "Cdict['C11'] = uc.set_in_units(7.436, '10^11 * dyn/cm^2')\n",
    "Cdict['C12'] = uc.set_in_units(2.619, '10^11 * dyn/cm^2')\n",
    "Cdict['C33'] = uc.set_in_units(5.160, '10^11 * dyn/cm^2')\n",
    "Cdict['C44'] = uc.set_in_units(3.135, '10^11 * dyn/cm^2')\n",
    "Cdict['C13'] = uc.set_in_units(2.917, '10^11 * dyn/cm^2')\n",
    "Cdict['C14'] = uc.set_in_units(1.541, '10^11 * dyn/cm^2')\n",
    "\n",
    "# Set values to existing C using the rhombohedral function\n",
    "C.rhombohedral(**Cdict)\n",
    "\n",
    "# Show that Cij array is properly constructed\n",
    "print('Cij (GPa) ->')\n",
    "print(uc.get_in_units(C.Cij, 'GPa'))\n",
    "print()\n",
    "\n",
    "# Show that 2 * C66 = C11 - C12\n",
    "print('Does 2*C66 = C11-C12?')\n",
    "print('C11 - C12 (GPa) =', uc.get_in_units(C.Cij[0,0] - C.Cij[0,1], 'GPa'))\n",
    "print('2 * C66 (GPa) =  ', uc.get_in_units(2 * C.Cij[5,5], 'GPa'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. ElasticConstants.tranform()<a id='section4'></a>\n",
    "\n",
    "The transform() method is included for convenience allowing for the elastic constants to be transformed to a different set of axes.\n",
    "\n",
    "Parameters\n",
    "\n",
    "- **axes** (3, 3) array giving three right-handed orthogonal vectors to use for transforming.\n",
    "- **tol** optional relative tolerance to use in identifying near-zero terms.\n",
    "\n",
    "Returns a new ElasticConstants object with constants transformed to the new axes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cij (GPa) ->\n",
      "[[176.2 124.9 124.9   0.    0.    0. ]\n",
      " [124.9 176.2 124.9   0.    0.    0. ]\n",
      " [124.9 124.9 176.2   0.    0.    0. ]\n",
      " [  0.    0.    0.   81.8   0.    0. ]\n",
      " [  0.    0.    0.    0.   81.8   0. ]\n",
      " [  0.    0.    0.    0.    0.   81.8]]\n",
      "\n",
      "After transforming, Cij (GPa) ->\n",
      "[[232.35  68.75 124.9    0.     0.     0.  ]\n",
      " [ 68.75 232.35 124.9    0.     0.     0.  ]\n",
      " [124.9  124.9  176.2    0.     0.     0.  ]\n",
      " [  0.     0.     0.    81.8    0.     0.  ]\n",
      " [  0.     0.     0.     0.    81.8    0.  ]\n",
      " [  0.     0.     0.     0.     0.    25.65]]\n"
     ]
    }
   ],
   "source": [
    "# Set C back to a cubic system (Cu this time)\n",
    "Cdict = {}\n",
    "Cdict['C11'] = uc.set_in_units(1.762, '10^12 * dyn/cm^2')\n",
    "Cdict['C12'] = uc.set_in_units(1.249, '10^12 * dyn/cm^2')\n",
    "Cdict['C44'] = uc.set_in_units(0.818, '10^12 * dyn/cm^2')\n",
    "\n",
    "C.cubic(**Cdict)\n",
    "print('Cij (GPa) ->')\n",
    "print(uc.get_in_units(C.Cij, 'GPa'))\n",
    "print()\n",
    "\n",
    "\n",
    "# Transform by a 45 degree rotation around z-axis\n",
    "axes = [[ 1, 1, 0],\n",
    "        [-1, 1, 0],\n",
    "        [ 0, 0, 1]]\n",
    "newC = C.transform(axes)\n",
    "\n",
    "print('After transforming, Cij (GPa) ->')\n",
    "print(uc.get_in_units(newC.Cij, 'GPa'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Shear and bulk modulus estimates<a id='section5'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.1. ElasticConstants.bulk()\n",
    "\n",
    "Three style options of estimates are available:\n",
    "\n",
    "- **'Voigt'** Voigt estimate.\n",
    "\n",
    "$$ K_{Voigt} = \\frac{ \\left(C_{11} + C_{22} + C_{33} \\right) + 2 \\left(C_{12} + C_{13} + C_{23} \\right) }{9}$$\n",
    "\n",
    "- **'Reuss'** Reuss estimate.\n",
    "\n",
    "$$ K_{Reuss} = \\frac{1}{ \\left( S_{11} + S_{22} + S_{33} \\right) + 2 \\left(S_{12} + S_{13} + S_{23} \\right) }$$\n",
    "\n",
    "- **'Hill'** Hill estimate (default).\n",
    "\n",
    "$$ K_{Hill} = \\frac{K_{Reuss} + K_{Voigt}}{2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Voigt bulk modulus estimate = 142.0 GPa\n",
      "Reuss bulk modulus estimate = 141.9999999999999 GPa\n",
      "Hill bulk modulus estimate =  141.99999999999994 GPa\n"
     ]
    }
   ],
   "source": [
    "print('Voigt bulk modulus estimate =', uc.get_in_units(C.bulk('Voigt'), 'GPa'), 'GPa')\n",
    "print('Reuss bulk modulus estimate =', uc.get_in_units(C.bulk('Reuss'), 'GPa'), 'GPa')\n",
    "print('Hill bulk modulus estimate = ', uc.get_in_units(C.bulk(), 'GPa'), 'GPa')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2. ElasticConstants.shear()\n",
    "\n",
    "Three style options of estimates are available:\n",
    "\n",
    "- **'Voigt'** Voigt estimate.\n",
    "\n",
    "$$ \\mu_{Voigt} = \\frac{ \\left( C_{11} + C_{22} + C_{33} \\right) - \\left( C_{12} + C_{23} + C_{13} \\right) + \n",
    "   3 \\left( C_{44} + C_{55} + C_{66} \\right) }{15} $$\n",
    "\n",
    "- **'Reuss'** Reuss estimate.\n",
    "\n",
    "$$ \\mu_{Reuss} = \\frac{15}{4 \\left( S_{11} + S_{22} + S_{33} \\right) - 4 \\left( S_{12} + S_{23} + S_{13} \\right) +\n",
    "  3 \\left( S_{44} + S_{55} + S_{66} \\right)} $$\n",
    "\n",
    "- **'Hill'** Hill estimate (default).\n",
    "\n",
    "$$ \\mu_{Hill} = \\frac{\\mu_{Reuss} + \\mu_{Voigt}}{2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Voigt shear modulus estimate = 59.34000000000001 GPa\n",
      "Reuss shear modulus estimate = 43.611930991477855 GPa\n",
      "Hill shear modulus estimate =  51.47596549573894 GPa\n"
     ]
    }
   ],
   "source": [
    "print('Voigt shear modulus estimate =', uc.get_in_units(C.shear('Voigt'), 'GPa'), 'GPa')\n",
    "print('Reuss shear modulus estimate =', uc.get_in_units(C.shear('Reuss'), 'GPa'), 'GPa')\n",
    "print('Hill shear modulus estimate = ', uc.get_in_units(C.shear(), 'GPa'), 'GPa')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Crystal system normalized values<a id='section6'></a>\n",
    "\n",
    "The elastic constants tensor should have certain components that are equal or dependent based on the crystal symmetry of the system.  However, calculations of the elastic constants for an atomic system may show some variability across the measured values of these symmetrically-dependent components.  The methods listed here help to handle the symmetry components."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.1. ElasticConstants.normalized_as()\n",
    "\n",
    "Returns a new ElasticConstants object where values of the current are averaged or zeroed out according to a standard crystal system setting.\n",
    "\n",
    "**NOTE:** no validation checks are made to evaluate whether such normalizations should be done! That is left up to you (compare values before and after normalization).\n",
    "\n",
    "Parameters\n",
    "\n",
    "- **crystal_system** (*str*) Indicates the crystal system representation to use when building a data model.\n",
    "\n",
    "Returns\n",
    "\n",
    "- (*atomman.ElasticConstants*) The elastic constants normalized according to the crystal system symmetries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C_with_v:\n",
      "[[176.1993 124.9016 124.9004   0.0036  -0.0042   0.003 ]\n",
      " [124.9016 176.1998 124.9007   0.0003   0.      -0.0041]\n",
      " [124.9004 124.9007 176.2035  -0.0011  -0.0004   0.0012]\n",
      " [  0.0036   0.0003  -0.0011  81.7983   0.0027   0.0006]\n",
      " [ -0.0042   0.      -0.0004   0.0027  81.7951   0.0015]\n",
      " [  0.003   -0.0041   0.0012   0.0006   0.0015  81.8033]]\n",
      "\n",
      "C_norm\n",
      "[[176.2009 124.9009 124.9009   0.       0.       0.    ]\n",
      " [124.9009 176.2009 124.9009   0.       0.       0.    ]\n",
      " [124.9009 124.9009 176.2009   0.       0.       0.    ]\n",
      " [  0.       0.       0.      81.7989   0.       0.    ]\n",
      " [  0.       0.       0.       0.      81.7989   0.    ]\n",
      " [  0.       0.       0.       0.       0.      81.7989]]\n"
     ]
    }
   ],
   "source": [
    "# Create a new ElasticConstants object with some variability\n",
    "scale = uc.set_in_units(1e-2, 'GPa')\n",
    "Cij_with_v = C.Cij + scale * np.random.rand(6,6) - scale/2\n",
    "for i in range(6):\n",
    "    for j in range(i, 6):\n",
    "        Cij_with_v[i,j] = Cij_with_v[j, i] = (Cij_with_v[i,j]+Cij_with_v[j, i])/2\n",
    "C_with_v = am.ElasticConstants(Cij=Cij_with_v)\n",
    "\n",
    "print('C_with_v:')\n",
    "print(uc.get_in_units(C_with_v.Cij, 'GPa'))\n",
    "print()\n",
    "\n",
    "# Normalize to cubic\n",
    "C_norm = C_with_v.normalized_as('cubic')\n",
    "print('C_norm')\n",
    "print(uc.get_in_units(C_norm.Cij, 'GPa'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.2. is_normal()\n",
    "\n",
    "Checks if current elastic constants agree with values normalized to a specified crystal family (within tolerances).\n",
    "        \n",
    "Parameters\n",
    "        \n",
    "- **crystal_system** (*str*) Indicates the crystal system representation to use when building a data model.\n",
    "- **atol** (*float, optional*) Absolute tolerance to use.  Default value is 1e-4.\n",
    "- **rtol** (*float, optional*) Relative tolerance to use.  Default value is 1e-4.\n",
    "        \n",
    "Returns\n",
    "\n",
    "- (*bool*) True if all Cij match within the tolerances, false otherwise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C_with_v.is_normal('cubic', rtol=1e-4, atol=1e-4) -> True\n",
      "C_with_v.is_normal('cubic', rtol=1e-7, atol=1e-7) -> False\n",
      "C_with_v.is_normal('hexagonal') -> False\n",
      "\n",
      "C_norm.is_normal('cubic', rtol=1e-4, atol=1e-4) -> True\n",
      "C_norm.is_normal('cubic', rtol=1e-7, atol=1e-7) -> True\n",
      "C_norm.is_normal('hexagonal') -> False\n"
     ]
    }
   ],
   "source": [
    "print(\"C_with_v.is_normal('cubic', rtol=1e-4, atol=1e-4) ->\", C_with_v.is_normal('cubic', rtol=1e-4, atol=1e-4))\n",
    "print(\"C_with_v.is_normal('cubic', rtol=1e-7, atol=1e-7) ->\", C_with_v.is_normal('cubic', rtol=1e-7, atol=1e-7))\n",
    "print(\"C_with_v.is_normal('hexagonal') ->\", C_with_v.is_normal('hexagonal'))\n",
    "print()\n",
    "\n",
    "print(\"C_norm.is_normal('cubic', rtol=1e-4, atol=1e-4) ->\", C_norm.is_normal('cubic', rtol=1e-4, atol=1e-4))\n",
    "print(\"C_norm.is_normal('cubic', rtol=1e-7, atol=1e-7) ->\", C_norm.is_normal('cubic', rtol=1e-7, atol=1e-7))\n",
    "print(\"C_norm.is_normal('hexagonal') ->\", C_norm.is_normal('hexagonal'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. ElasticConstants.model()<a id='section7'></a>\n",
    "\n",
    "The elastic constants can also be saved/retrieved as a JSON/XML data model. This is useful as it captures not only the elastic constant values but also the associated units.\n",
    "\n",
    "Parameters\n",
    "\n",
    "- **model** (*DataModelDict, string, or file-like object, optional*) Data model containing exactly one 'elastic-constants' branch to read.\n",
    "- **unit** (*str, optional*) Units or pressure to save values in when building a data model. Default value is None (no conversion).\n",
    "- **crystal_system** (*str, optional*) Indicates the crystal system representation to use when building a data model. Default value is 'triclinic' (save all values in Cij).\n",
    "        \n",
    "Returns \n",
    "\n",
    "- (*DataModelDict*) If model is not given as a parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{\n",
      "    \"elastic-constants\": {\n",
      "        \"Cij\": {\n",
      "            \"value\": [\n",
      "                176.20000000000002,\n",
      "                124.9,\n",
      "                124.9,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                124.9,\n",
      "                176.20000000000002,\n",
      "                124.9,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                124.9,\n",
      "                124.9,\n",
      "                176.20000000000002,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                81.8,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                81.8,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                0.0,\n",
      "                81.8\n",
      "            ],\n",
      "            \"shape\": [\n",
      "                6,\n",
      "                6\n",
      "            ],\n",
      "            \"unit\": \"GPa\"\n",
      "        }\n",
      "    }\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "# Generate data model based on C\n",
    "model = C.model(unit='GPa', crystal_system='cubic')\n",
    "\n",
    "# Show json version of model\n",
    "print(model.json(indent=4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{\"elastic-constants\": {\"Cij\": {\"value\": [176.20000000000002, 124.9, 124.9, 0.0, 0.0, 0.0, 124.9, 176.20000000000002, 124.9, 0.0, 0.0, 0.0, 124.9, 124.9, 176.20000000000002, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 81.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 81.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 81.8], \"shape\": [6, 6], \"unit\": \"GPa\"}}}\n"
     ]
    }
   ],
   "source": [
    "# If crystal_system is not given, model will contain all 21 triclinic constants\n",
    "model = C.model(unit='GPa')\n",
    "print(model.json())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model can also be read in to an ElasticConstants object either during initialization or using the model() method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[176.2 124.9 124.9   0.    0.    0. ]\n",
      " [124.9 176.2 124.9   0.    0.    0. ]\n",
      " [124.9 124.9 176.2   0.    0.    0. ]\n",
      " [  0.    0.    0.   81.8   0.    0. ]\n",
      " [  0.    0.    0.    0.   81.8   0. ]\n",
      " [  0.    0.    0.    0.    0.   81.8]]\n"
     ]
    }
   ],
   "source": [
    "C = am.ElasticConstants(model=model)\n",
    "print(uc.get_in_units(C.Cij, 'GPa'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[176.2 124.9 124.9   0.    0.    0. ]\n",
      " [124.9 176.2 124.9   0.    0.    0. ]\n",
      " [124.9 124.9 176.2   0.    0.    0. ]\n",
      " [  0.    0.    0.   81.8   0.    0. ]\n",
      " [  0.    0.    0.    0.   81.8   0. ]\n",
      " [  0.    0.    0.    0.    0.   81.8]]\n"
     ]
    }
   ],
   "source": [
    "newC.model(model=model)\n",
    "print(uc.get_in_units(newC.Cij, 'GPa'))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
